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Abstract 


A  system  of  ordinary  differential  equations  is  formulated  to  describe  the  pathogenesis  of  HIV 
infection,  wherein  certain  important  features  that  have  been  shown  important  by  recent  ex¬ 
perimental  research  are  incorporated  in  the  model.  These  include  the  role  of  CD4+  memory 
cells  that  serve  as  a  major  reservoir  of  latently  infected  cells,  a  critical  role  for  T-helper  cells 
in  the  generation  of  CD8  memory  cells  capable  of  efficient  recall  response,  and  stimulation 
by  antigens  other  than  HIV.  A  stability  analysis  illustrates  the  capability  of  this  model  in 
admitting  multiple  locally  asymptotically  stable  (locally  a.s.)  off-treatment  equilibria.  The 
phenomenon  of  “viral  blips”  experienced  by  some  patients  on  therapy  with  viral  load  levels 
suppressed  below  the  detection  limit  is  also  investigated.  Censored  clinical  data  is  used  to 
demonstrate  that  this  model  provides  reasonable  fits  to  all  the  patient  data  available  for  this 
study  and,  moreover,  that  it  exhibits  impressive  predictive  capability. 
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1  Introduction 


Since  the  seminal  work  of  Ho,  et  ah,  [18]  demonstrated  the  promise  for  elucidating  HIV 
disease  mechanisms  through  mathematical  modeling,  a  wide  variety  of  models  have  been 
proposed  to  describe  various  aspects  of  in- host  HIV  infection  dynamics  (e.g.,  [1,  2,  3,  12,  16, 
25,  26]).  The  most  basic  of  these  models  typically  include  two  or  three  of  the  key  dynamic 
compartments:  virus,  uninfected  target  cells,  and  infected  cells.  These  compartmental  de¬ 
pictions  lead  to  systems  of  linear  or  nonlinear  ordinary  differential  equations  in  terms  of  state 
variables  representing  the  concentrations  in  each  compartment  and  parameters  describing  vi¬ 
ral  production  and  clearance,  cell  infection  and  death  rate,  treatment  efficacy,  etc.  Solutions 
for  the  model  states  yield  the  time  course  of  viral  load  and  CD4+  counts,  for  example. 

Although  such  models  can  be  expected  only  to  approximate  the  myriad  processes  underlying 
HIV  pathogenesis,  when  used  in  conjunction  with  data  as  part  of  designed  experiments,  these 
models  can  be  powerful  tools  in  answering  questions  about  the  pathogenesis  of  HIV  infection 
or  similar  biological  processes.  Mathematical  models  can  also  stimulate  further  clinical  and 
laboratory  research  [26].  For  example,  early  applications  of  linear  systems  to  short-term 
data  on  patients  undergoing  ARV  therapy  suggested  the  now  widely-held  theory  of  very 
rapid  and  constant  turnover  of  viral  and  infected  cell  populations  [18,  25],  contradicting 
previous  assumptions  that  stable  viral  and  CD4+  concentrations  during  the  clinical  latency 
period  of  chronic  HIV  infection  are  clue  to  absence  of  significant  viral  replication. 

Some  important  features  of  HIV  pathogenesis  and  the  cellular  immune  response  that  have 
emerged  in  recent  research  include  a  clearer  delineation  of  the  importance  of  memory  CD4+ 
T-cells  as  a  latent  reservoir  of  HIV  [17,  33]  and  a  critical  role  for  T-helper  cells  in  the 
generation  of  CD8  memory  cells  capable  of  an  efficient  recall  response  [4,  8,  19].  The  authors 
in  [33]  indicate  that,  even  in  treated  patients  who  have  had  no  detectable  viremia  for  as  long 
as  7  years,  the  reservoir  decays  so  slowly  that  early  initiation  of  Highly  Active  Anti- Retroviral 
Therapy  (HAART)  with  the  goal  of  virus  eradication  is  not  likely  to  succeed.  This  motivates 
us  to  develop  a  model  that  incorporates  these  features.  In  any  discussions  of  mathematical 
modeling  of  complex  systems  it  is  appropriate  to  point  out  that  while  complex  models  may 
be  needed  to  provide  accurate  descriptions  of  the  underlying  dynamics,  the  models  are  most 
useful  when  they  can  be  compared  to  clinical  and/or  experimental  data  and  can  also  be 
used  for  prediction.  In  developing  models  for  HIV  infection  and  treatment  or  some  other 
biological  phenomenon,  this  requires  a  balance  between  complexity  and  utility.  Hence,  in 
this  paper  we  do  not  try  to  formulate  a  model  that  reflects  all  features  of  cellular  immune 
response  as  well  as  all  host  and  viral  factors.  Instead,  we  attempt  to  develop  a  model  that 
can  capture  the  most  salient  features  of  disease  progression,  that  can  be  supported  and 
validated  based  on  data,  one  for  which  parameters  can  be  plausibly  estimated,  one  that  has 
predictive  capabilities,  and  one  for  which  control/drug  therapy  design  is  tractable.  While  the 
model  developed  and  analyzed  here  is  new,  it  modifies  and  extends  both  conceptually  and 
structurally  the  predictive  model  in  [3].  That  model  included  both  CD4+/viral  dynamics 
as  in  models  discussed  in  [12]  as  well  as  immune  response  compartments  whose  importance 
have  been  earlier  established  [10,  24,  39]  -  see  the  discussions  in  [1], 
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The  paper  is  organized  as  follows.  In  Section  2,  a  system  of  ordinary  differential  equations  is 
developed  to  describe  the  pathogenesis  of  HIV  and  the  cellular  immune  response.  In  Section 
3  we  discuss  the  ability  of  the  model  to  admit  multiple  locally  a.s.  off-treatment  equilibria. 
In  Section  4  we  investigate  a  possible  mechanism  for  producing  viral  blips  and  elucidate  the 
role  of  latently-infected  CD4+  memory  cells  and  the  effect  of  CD4+  help  on  CD8+  memory 
during  the  ensuing  immune  response.  The  expectation  maximization  algorithm  leading  to 
weighted  least-squares  techniques  is  employed  in  Section  5  to  fit  the  model  to  clinical  data 
with  lower  limit  censoring.  The  predictive  capability  of  the  model  is  also  investigated  by  using 
simulation  results,  with  parameters  estimated  from  only  half  of  the  longitudinal  observations, 
to  predict  the  immune  response  in  the  latter  half  and  comparing  these  predictions  to  clinical 
observations  and  model  results  obtained  from  fitting  the  full  longitudinal  data  sets.  Finally 
we  close  with  summary  conclusions  and  remarks  in  Section  6. 


2  HIV  Model 


The  model  we  develop  in  this  paper  conceptually  modifies  and  extends  the  model  in  [3], 
wherein  two  types  of  target  cells  (CD4+  T-cclls  and  macrophages),  along  with  their  cor¬ 
responding  infected  states,  free  virus,  and  immune  effector  cells  (CTL)  are  included  in  the 
model.  Clinical  data  fitting  results  show  that  the  preliminary  model  of  [3]  provides  rea¬ 
sonable  fits  to  most  patient  data  and  has  impressive  predictive  capability  when  comparing 
model  simulations,  with  parameters  based  on  estimation  using  only  half  of  the  longitudinal 
observations,  to  the  full  longitudinal  data  sets.  However,  that  model  does  not  incorporate 
some  important  features  of  HIV  pathogenesis  and  the  cellular  immune  response,  such  as 
CD4+  memory  cells  as  the  major  reservoir  of  latently  infected  cells  and  a  critical  role  for 
T-helper  cells  in  the  generation  of  CD8  memory  cells  capable  of  an  efficient  recall  response. 
To  incorporate  these  important  features,  we  thus  seek  a  model  that  includes  some  measure 
of  CD4+  T-helper  cells,  infected  memory  CD4+  T-cells  and  HIV-specific  memory  CD8+ 
T-cells.  To  retain  the  simplicity  of  the  model,  secondary  target  cells,  such  as  macrophages, 
are  not  included  as  a  compartment  in  our  new  model.  It  is  worth  noting  that  omitting  the 
secondary  target  cells  should  not  affect  our  clinical  data  fitting  and  predictive  capabilities 
since  this  type  of  cell,  even  though  it  is  very  important  at  the  beginning  of  infection,  does 
not  contribute  significantly  to  the  virus  pool  in  the  long  run.  The  model  compartments  are 
illustrated  in  Table  1,  wherein  the  resting  CD4+  T-cells  (To)  are  assumed  to  include  naive 
CD4+  T-cells  and  memory  CD4+  T-cclls.  This  is  reasonable  since  these  two  types  of  cells 
have  similar  behavior  such  as  longer  life  spans  and  distribution  in  the  lymphoid  tissue.  Once 
these  resting  CD4+  T-cells  become  activated,  either  through  antigen  priming  of  naive  cells 
or  reactivation  of  memory  cells,  they  are  more  susceptible  to  HIV  infection  than  resting  cells 
and  suffer  elevated  mortality.  Hence,  we  include  these  activated  naive  cells  and  reactivated 
memory  cells  in  the  other  compartment  as  activated  CD4+  T-cells  (Tf).  Infected  resting 
and  activated  cells  are  represented  by  the  T2*  and  7j*  states,  respectively.  A  schematic  of  this 
new  model  is  depicted  in  Fig.  1. 

The  corresponding  compartmental  ordinary  differential  equation  (ODE)  model  for  in-host 
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states 

unit 

description 

T] 

cells//il-blood 

uninfected  activated  CD4+  T-cclls 

T* 

cells//il-blood 

infected  activated  CD4+  T-cells 

t2 

cells//il-blood 

uninfected  resting  CD4+  T-cells 

t; 

cells/ /il-blood 

infected  resting  CD4+  T-cells 

Vj 

RNA  copies/ml-plasma 

infectious  free  virus 

Vni 

RNA  copies/ml-plasma 

non-infections  free  virus 

Ei 

cells/ /il-blood 

HIV-specific  effector  CD8+  T-cclls 

E2 

cells/ /il-blood 

HIV-specific  memory  CD8+  T-cells 

Tabic  1:  Model  States. 


Figure  1:  Flow  chart  of  model  (2.1),  where  PI  and  RTI  denote  protease  inhibitor  and  reverse 
transcriptase  inhibitor,  respectively. 
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HIV  infection  dynamics  is  based  on  balance  laws  and  is  given  by 


T\  —  —d\Ti  —  (1  —  £i(t))kiViTi  —  7 tTi  +  Pt  ^ vi'+kv  ^  °A)  ^2’ 

T*  =  (1  -  UmiVjTi  -  ST *  -  mE{T{  -  7tT*  +  pT  T * 

r2  =  +  1ttx  -  d2r2  -  (1  -  +  aA)  t2, 

t;  =  7 tT*  +  (1  -  Mi(t))k2V tT2  -  d2T*  -  (V;^  +  aA)  T*, 

W  =  (1  -  6W)103^Vt^  -  cV7  -  103[(1  -  6(*))PiM  +  (1  -  W)P2^2T2]VJ, 


Viv/  =  Ut)W3NTST*  -  cVni , 


-Ei  —  A^  + 


fesiT* 

E+Kbl 


Ei 


dET* 

T*+Kd 


E\  —  8  ei  Ei  —  7  e 


Ti+T* 

n+T*+K. 


Ei  +  pgagVf  Eo 
V,+Kv  ^2, 


-E2  —  7e 


ri  +'ri*  p  I 

Ti+Tf+AT^  1  ' 


bE2Kb2 

E2+Kb2 


E2  —  8eiE2 


aEVi 

V/+AV 


E 


2, 


(2.1) 


with  an  initial  condition  vector 


[Ti(0),  T*(0),T2(0),  T*( 0),  V,(0),  Vv/(0),  -Ei(O),  £2(0)]T. 

Here  the  factors  103  are  introduced  to  convert  between  microliter  and  milliliter  scales,  pre¬ 
serving  the  units  from  some  of  the  earlier  published  papers  [1,  10].  The  treatment  factors 
£i (t)  =  ei u(t)  and  £2(£)  =  e2u(t)  represent  the  effective  treatment  impact,  consisting  of  effi¬ 
cacy  factors  ei  modeling  the  relative  effectiveness  of  reverse  transcriptase  inhibitor  (RTI),  e2 
describing  the  relative  effectiveness  of  protease  inhibitor  (PI),  and  a  time-dependent  treat¬ 
ment  function  u(t)  (0  <  u(t)  <  1)  representing  HAART  drug  level,  where  u(t)  =  0  is  fully  off 
and  u(t)  =  1  is  fully  on.  Since  HIV  treatment  is  nearly  always  administered  as  combination 
therapy,  we  do  not  consider  the  possibility  of  monotherapy,  even  for  a  limited  period  of  time, 
though  this  could  be  implemented  by  considering  separate  treatment  functions. 

The  input  term  A t  ^  l\sK  for  the  T2  compartment  is  used  to  account  for  the  source  rate  of 
uninfected  resting  CD4+  T-cells.  This  term  depends  on  the  viral  load  level  since  the  thymus 
production  can  be  diminished  if  the  viral  load  is  too  high  [22],  To  limit  the  introduction 
of  additional  parameters,  we  assume  that  uninfected  and  infected  resting  CD4+  T-cells  (T2 
and  T-2,  respectively)  have  the  same  natural  death  rate  d2.  We  remark  that  activated  CD4+ 
T-cells  have  a  higher  natural  death  rate  than  resting  memory  and  naive  cells,  and  we  use 
d\  to  denote  the  natural  death  rate  of  uninfected  activated  CD4+  T-cells  T\.  The  immune 
effector  cells  E\  remove  infected  activated  cells  CD4+  T-cells  T*  from  the  system  by  cell  lysis 
with  a  rate  m.  However,  immune  effector  cells  do  not  remove  infected  resting  cells  ,  since 
these  cells  are  in  a  quiescent  state  where  the  virus  is  not  replicating  and,  thereby,  escape  the 
detection  of  the  immune  effector  cells.  These  infected  resting  cells  are  assumed  to  become 
targets  for  lysis  only  after  activation  [9]. 

The  infected  activated  cells  T1*  result  from  encounters  between  uninfected  activated  cells  T\ 
and  free  infectious  virus  Vj  with  infection  rate  k\ .  The  resulting  term  k\V\T\  is  modified  by  a 
factor  1  —  £i(t)  to  account  for  RTI  treatment.  Infection  of  the  resting  T-ccll  compartment  T2, 
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which  is  comprised  of  both  memory  and  naive  CD4+  T-cells,  can  occur  in  a  number  of  ways. 
First,  the  most  commonly  transmitted  R5  virus  form  of  HIV-1  that  utilize  the  chemokine 
receptor  CCR5  can  enter  a  subset  of  resting  memory  cells  that  express  sufficient  levels  of 
CCR5  to  support  infection  [9].  In  addition,  the  X4  form  of  the  virus  can  infect  resting  CD4+ 
T-cells,  whether  they  belong  to  the  naive  or  memory  subsets.  However,  infection  of  naive 
and  memory  cells  through  these  routes  occurs  much  less  frequently  than  infection  of  Ti,  and, 
once  infected,  these  cells  often  do  not  progress  to  a  long-term  stably-infected  state  in  which 
the  virus  is  integrated  into  the  host  DNA.  In  addition,  it  has  been  shown  that  infected  naive 
CD4+  T-cells  do  not  significantly  contribute  to  the  pool  of  infected  resting  CD4+  T-cells 
[11].  Hence,  the  term  (1  —  /£i(t))&2V/T2  is  used  to  represent  the  infection  process  that  results 
from  encounters  between  the  uninfected  resting  CD4+  T-cells  and  free  virus  Vf,  but  with  an 
infection  rate  fc2  <  k\  to  account  for  a  significantly  lower  rate  of  infection  as  compared  to 
activated  CD4+  T-cells.  The  treatment  factor  £i(t)  is  potentially  more  effective  in  Ti  than 
in  T2,  where  the  efficacy  is  modelled  by  f£i(t)  with  0  <  /  <  1. 

A  much  more  stable  form  of  latent  infection  arises  when  activated  CD4+  T-cells  that  have 
integrated  HIV-1  DNA  survive  long  enough  to  revert  back  to  resting  memory  state,  and 
latently  infected  resting  CD4+  T-cells  with  integrated  HIV-1  DNA  are  present  in  all  infected 
individuals  but  only  at  low  frequency  [9,  14].  Hence,  the  terms  involving  7 pT i*  are  included 
in  the  model  to  account  for  the  phenomenon  of  differentiation  of  infected  activated  CD4+ 
T-cells  into  infected  memory  or  resting  CD4+  T-cells  T2*  at  rate  7 •?.  For  simplicity,  the  rate 
at  which  uninfected  activated  CD4+  T-cells  T\  differentiate  into  uninfected  resting  CD4+ 
T-cells  T-2  is  also  assumed  to  be  7^;  the  model  could  be  extended  easily  to  the  case  with 
different  differentiation  rates. 

As  the  authors  of  [9]  concluded,  there  is  turnover  in  the  latent  reservoir  when  patients  are 
viremic,  but  the  degree  of  turnover  depends  on  the  level  of  viremia.  We  thus  assume  that 
the  activation  of  infected  HIV-specific  resting  CD4+  T-cells  T2*  depends  on  the  virus  con¬ 
centration  with  a  half-saturation  constant  Kv.  Hence,  the  terms  involving  T2*  are 

used  to  represent  the  activation  of  infected  HIV-specific  resting  CD4+  T-cells  with  maxi¬ 
mum  activation  rate  of  «t-  Again  to  preserve  the  simplicity  of  this  model,  we  assume  that 
activation  of  uninfected  HIV-specific  resting  CD4+  T-cells  T2  also  depends  on  the  virus  con¬ 
centration,  with  a  half  saturation  constant  KVl  and  that  the  maximum  activation  rate  is  also 
aT .  Thus,  the  terms  involving  yI^v  T2  represent  the  activation  of  uninfected  HIV-specific 
resting  CD4+  T-cells.  In  order  to  incorporate  the  activation  of  resting  CD4+  T-cells  by 
some  non-HIV  antigen  and  preserve  the  simplicity  of  the  model,  we  include  the  simple  terms 
oaT2  and  a^T2  into  our  model  to  describe  this  phenomenon,  with  cla  being  the  activation 
rate  by  non-HIV  antigen.  The  parameter  here  can  be  utilized  as  a  constant  to  represent 
a  chronic  level  of  infection  or  as  a  function  of  time  t  to  describe  infections  that  are  cleared 
by  the  body.  These  activation  terms  represent  losses  to  the  T2  and  T2  compartments,  with 
corresponding  gain  terms  for  the  T\  and  T*  compartments.  However,  the  gain  terms  for  T\ 
and  T*  include  a  multiplicative  factor  pt  to  account  for  the  net  proliferation  due  to  clonal 
expansion  and  programmed  contraction.  For  simplicity,  we  assume  that  uninfected  and  in¬ 
fected  CD4+  T-cells  have  the  same  expansion  factor  p?',  again  this  can  be  readily  extended 
to  include  processes  with  different  expansion  factors. 
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Virus  in  the  reservoir  T2*  of  infected  resting  CD4+  T-cells  is  latent  and  no  virus  can  be 
produced  by  these  cells  unless  they  are  activated  [9].  Hence,  free  virus  particles  Vj  are 
produced  only  by  activated  infected  CD4+  T-cells  during  viral  budding  leading  up  to  viral 
produced  lysis  5T±  of  the  CD4+  T-cells.  The  parameter  NT  accounts  for  the  number  of 
RNA  copies  produced  during  this  process  in  the  viral  source  term  (1  —  £2(t))103 NE5Tl . 
In  addition  to  a  natural  clearance  rate  c,  we  also  include  terms  103[(1  —  £i(£)) pikiTi  and 
(1  —  f  £,i(t))  p2k2T2]Vj  in  the  free  virus  compartment  Vj  to  account  for  the  removal  of  free  virus 
that  takes  place  when  free  virus  infects  7\  and  T2.  We  make  the  simplifying  assumption  that 

Pi  =  1  C°^  6S  — - — ,  i.e.,  one  free  virus  particle  is  responsible  for  each  new  infection.  This 
cells  ml-plasma 

could  be  adapted  easily  for  multiple  virus  particles  being  responsible  for  each  new  infection 
by  choosing  p,  >  1.  Since  clinical  measurements  of  viral  load  do  not  differentiate  between 
infectious  and  non-infectious  virus,  we  include  a  compartment  in  the  model  for  tracking  the 
amount  of  non-infectious  virus  Vni ■  The  action  of  a  protease  inhibitor,  resulting  in  the 
production  of  non-infectious  virus  Vni  by  infected  cells  is  modeled  by  £2.  It  should  be  noted 
that  the  inclusion  of  this  additional  state  does  not  affect  the  dynamics  of  the  other  state 
variables. 


The  source  term  XE,  the  constant  death  term  SE i,  and  the  nonlinear  infected  cell-dependent 


birth 


Tf+Kbl 


Ei  and  death  E\  terms  in  the  E\  compartment  are  adopted  from  the 


model  in  [2,  3],  where  the  authors  suggested  that,  by  including  such  terms  in  the  immune 
effector  compartment,  the  model  can  admit  multiple  stable  off-treatment  steady  states  and 
exhibit  transfer  between  “healthy”  and  “unhealthy”  locally  stable  steady  states  via  optimal 
or  suboptimal  structure  treatment  interruptions  (STI)  therapies.  This  makes  it  a  good 
candidate  for  our  investigation.  Memory  CD8+  T-cells  are  also  subject  to  strict  homeostatic 
control  [34];  background  expansion  of  memory  cells  through  intermittent  cell  division  being 
countered  by  an  equivalent  level  of  cell  death.  Hence,  we  include  the  term  E2  —  5E2E2 

for  homeostatic  regulation  in  the  E2  compartment,  similar  to  that  used  in  [37].  In  the 
homeostatic  regulation,  bE2  represents  the  maximum  proliferation  rate  and  SE2  corresponds 
to  the  death  rate,  where  the  proliferation  signal  decreases  linearly  with  population  size. 


The  term  1et^t*+k  in  the  model  is  used  to  include  the  essential  role  that  activated 
CD4+  T-cells  play  in  the  generation  of  memory  CD8+  T-cells  during  the  priming  phase, 
where  parameter  K1  is  a  half-saturation  constant  and  ryE  is  the  maximum  rate  at  which 
Ei  differentiates  into  E-2.  Since  depletion  of  CD4+  cells  has  a  minimal  effect  during  the 
recall  response  [32,  38],  the  term  y^E2  for  reactivation  of  memory  CD8+  T-cells  is 
independent  of  CD4+  T-ccll  help.  Similar  to  the  activation  of  HIV-specific  resting  CD4+ 
T-cells,  we  assume  that  activation  of  HIV-specific  memory  CD8+  T-cells  also  depends  on 
the  virus  concentration.  For  simplicity,  we  use  the  same  half-saturation  constant  Kv  for 
the  activation  of  memory  CD8+  T-cells.  Since  CD8+  T-cells  tend  to  divide  sooner  and  to 
have  a  faster  rate  of  cell  division  than  CD4+  T-cells  [31],  we  use  a  different  parameter  pE 
to  account  for  the  net  proliferation  due  to  clonal  expansion  and  programmed  contraction  of 
activated  CD8+  T-cells  in  the  Ei  compartment. 
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3  Off-treatment  Stability  Analysis 


Our  model  choice  is  partly  motivated  by  its  admission  of  multiple  locally  a.s.  off-treatment 
equilibria  (corresponding  to  e\  =  0  and  €2  =  0),  which  is  important  in  the  sense  that  we  seek 
to  investigate  the  possibility  of  improvement  of  long-term  health  by  using  STI  during  acute 
infection  to  affect  a  change  from  an  “unhealthy”  equilibrium  to  a  “healthy”  one.  A  general 
analysis  of  the  equilibria  of  the  model  (2.1)  and  their  stability  properties  is  challenging 
due  to  the  complexity  of  the  system  and  is  generally  not  solvable  in  analytic  form.  Hence, 
in  this  section  we  illustrate  the  existence  of  multiple  locally  a.s.  equilibria  through  specific 
examples  and  numerically  investigate  the  behavior  of  these  equilibria  with  respect  to  changes 
in  parameters  and  initial  conditions. 

Uninfected  equilibrium  in  a  healthy  person  (a  a  =  0).  In  the  absence  of  non-HIV 

antigen  ( ctA  =  0),  it  is  easy  to  see  that  model  (2.1)  has  the  following  off-treatment  uninfected 
equilibrium: 

(0,  0,  A t/ d'2j  0,  0,  0,  A e/ Sei,  0).  (3.1) 

Substituting  the  above  equilibrium  into  the  Jacobian  matrix  of  model  (2.1),  we  find  that  the 
eigenvalues  of  this  matrix  are  given  by  —  (<5<5ei  +  ttiXe  +  7t£ei) / Sei,  — c,  —d.2,  —  g?2,  —  (cg?2  + 
IOOOA^AtVg^,  —  Sei,  i>E2  —  &E2  and  —d\  —  7 t-  Hence,  if  2  <  Se2  then  all  the  eigenvalues 
of  this  Jacobian  matrix  are  negative,  which  implies  that  under  this  case  equilibrium  (3.1)  is 
locally  a.s.  On  the  other  hand,  if  Se2  >  Se2,  then  we  obtain  a  different  uninfected  equilibrium: 

(0, 0,  A t/ d2,  0,  0, 0,  Xe/ Sei,  (&E2/ $E2  ~  1)  AJ^)-  (3.2) 

In  this  case,  we  find  that  the  eigenvalues  of  Jacobian  matrix  are  given  by  — (JJ^i  +  nrA e  + 
XtSei)  /  Sei,  —  c,  —di—'fTi  — (cg?2  +  1000/^2 At)/g?2,  —Sei,  ~ ^2(^2  —  Se2)/Se2,  —d>2  and  — c?2 - 
Hence,  equilibrium  (3.2)  is  also  locally  a.s.  Therefore,  for  the  case  <74  =  0,  model  (2.1) 
always  has  a  locally  a.s.  off-treatment  uninfected  equilibrium.  The  existence  of  a  locally 
a.s.  uninfected  equilibrium  is  biologically  reasonable  in  light  of  research  documenting  the 
existence  of  some  sex  workers  who,  though  regularly  exposed  to  HIV-contaminated  body 
fluids,  remain  HIV-negative  [28,  29].  We  note  that  this  feature  of  our  new  model  is  not 
present  in  the  earlier  model  investigated  in  [3]. 

Multiple  equilibria  in  a  healthy  person  [a a  =  0).  We  next  illustrate  the  existence 
of  multiple  equilibria  with  specific  examples  when  ola  =  0  and  the  values  of  all  the  other 
parameters  as  specified  in  Table  2,  unless  otherwise  stated.  Since  values  for  most  of  the 
parameters  in  this  model  can  not  be  found  in  the  literature,  the  values  listed  in  Table  2 
are  chosen  for  model  illustration  purposes.  When  determining  equilibria  for  this  complex 
model  there  are  usually  many  unstable  or  physically-meaningless  (eg.,  negative  state  values) 
equilibria.  In  the  discussion  that  follows,  unless  otherwise  stated,  we  focus  our  attention  on 
the  locally  a.s.,  physically-meaningful  equilibria  only. 

In  particular,  we  wish  to  explore  the  effect  that  CD4+  memory  cell  activation  plays  on 
the  equilibrium  states  available  to  the  system.  Since,  for  economy,  we  have  included  both 
naive  and  memory  cells  in  a  single  resting  CD4+  T-cell  compartment  T2  (or  T2*  for  the 


parameter  value 

parameter  value 

parameter  value 

d\  0.02  — — - 

day 

d  G  [0, 1] 

5  mi-piasnra 

10  .  _ 
copies  •  day 

5  0.7  — — - 

day 

nl-  blood 

m  0.01  — - - - — 

cells  •  day 

ceils 

ul-blood  •  day 

1 

d2  0.005  — 

day 

_  copies 

Ay  100 

ml-plasma 

AT  105 

ml-plasma 

/  0.34 

i 

7 t  0.005  - — 

day 

„„  copies  •  ml- blood 
Nt  100  cells  ■  ml-plasma 

e2  €  [0, 1] 

1 

c  13  — 

day 

cells 

Ae  °  001  ul-blood  ■  day 

bE\  0.3  - — 

day 

dE  0.25 

day 

ceils 

61  -1  ul-blood 

aA  [0' 11  ^ 

5  El  0.1- — 

day 

ceils 

d  °'5  ul-blood 

„  i 

aE  0.008  — — 

day 

Pt  1-2 

,  1  „  _9  mi-plasma 

k2  10  - T - y— 

copies  •  day 

1 

aE  0.1- - 

day 

PE  3 

ceils 

7  10  ul-blood 

7s  0.01  - — 

day 

bE2  0.001  — — 

day 

_ _ ceils 

Xb2  °°  ul-blood 

5e2  0.005  / 

day 

Table  2:  Parameter  values  used  in  model  (2.1)  to  illustrate  the  existence  of  multiple  equilib¬ 
ria.  Note  that  for  this  parameter  set  bE2  <  SE 2- 


state 

EQi 

Cirp 

8  x  10-3 

8  x  10-2 

1  x  10-3 

eq2 

eq3 

eq2 

eq3 

eq2 

T\ 

0 

184.3 

13.13 

352.1 

13.12 

40.27 

T* 

0 

0.05621 

5.330 

0.05582 

7.659 

0.05935 

t2 

1400 

731.2 

424.4 

142.8 

59.36 

1247 

rp* 

1  2 

0 

0.04393 

2.984 

0.005110 

0.4801 

0.1187 

Vj 

0 

265.1 

28410 

236.5 

40830 

310.0 

Vni 

0 

0 

0 

0 

0 

0 

E\ 

0.01 

799.2 

0.04574 

1422 

0.03872 

140.0 

e2 

0 

98.32 

0.002863 

184.5 

0.002520 

14.07 

Table  3:  Off-treatment  steady  states  for  model  (2.1)  with  =  0  day-1,  aE  =  8  x  10-3, 
8  x  10-2,  and  1  x  10-3  day-1,  and  the  values  of  all  the  other  parameters  as  specified  in  Table 
2.  Non-physical  and  unstable  steady  states  are  omitted. 
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infected  case),  the  activation  factor  represents  both  activation  of  naive  cells  and  the 

reactivation  of  memory  cells.  In  the  limit  of  large  Vj  the  maximum  activation  rate  is  ax- 
In  the  treatment  that  follows,  we  examine  the  role  of  ax  on  the  equilibrium  states  available 
to  the  system.  We  first  examine  the  off-treatment  equilibria  (listed  in  Table  3)  for  the  case 
ax  =  8  x  1CT3  day-1.  We  can  see  that  in  this  and  each  of  the  other  two  cases  we  consider, 
there  is  an  uninfected  equilibrium  similar  to  the  one  in  (3.1),  which  we  designate  as  EQ\. 
In  addition  to  the  uninfected  equilibrium,  we  have  two  infected  equilibria,  designated  EQ2 
and  EQ 3,  where  EQ2  represents  a  “healthy”  steady  state  with  immune  control  of  the  viral 
infection,  restoration  of  CD4+  T-cell  counts  (915.6  cells/pl-blood),  and  a  strong  CD8+ 
immune  response  (897.5  cells / /xl-blood) .  Equilibrium  EQ3  represents  an  “unhealthy”  steady 
state  corresponding  to  a  dangerously  high  viral  load  set  point,  lower  CD4+  T-cell  counts 
(445.8  cells//il-blood),  and  a  much  lower  CD8+  immune  response  (0.04860  cells/ //1-blood). 

If  we  begin  simulations  with  the  initial  conditions  (0,  0, 1400,  0,  z/,  0,  0.01,  0),  i.e.,  a  nontrival 
amount  v  of  infectious  free  virus,  then  the  solution  of  model  (2.1)  converges  to  either  the 
uninfected  equilibrium  EQi  or  the  “unhealthy”  infected  equilibrium  EQ3  (see  Fig.  2). 
The  particular  steady  state  that  the  model  converges  to  depends  on  the  viral  load  level  u. 
Simulation  results  reveal  that  the  solution  will  not  converge  to  EQ3  until  the  value  of  v  is 
close  to  6576  copies/ml-plasma. 


Figure  2:  Phase  diagram  showing  equilibrium  attained  as  a  function  of  the  initial  viral  load 
v  and  the  parameter  ax-  It  should  be  emphasized  that  this  plot  is  only  applicable  for  the 
particular  initial  condition  (0,  0, 1400,  0,  z/,  0,  0.01,  0). 

If  we  set  the  value  of  ax  larger,  such  as  8  x  10-2  day-1,  and  keep  the  values  of  all  the 
other  parameters  the  same,  then  we  still  have  three  locally  a.s.  off-treatment  equilibria: 
uninfected,  “healthy”  infected  and  “unhealthy”  infected,  where  the  uninfected  equilibrium 
is  the  same  as  before  (Table  3).  However,  the  “unhealthy”  infected  equilibrium  (EQ3)  in  this 
case  has  a  much  higher  viral  load,  a  much  lower  CD4+  T-cell  count  (80.62  cells///l-blood), 
and  a  lower  immune  response  (0.04124  cells/yul-blood)  than  the  EQ3  when  ax  =  8  x  10-3 
day-1.  Even  though  the  “healthy”  infected  equilibrium  EQ2  has  a  slightly  lower  viral  load 
level  and  a  much  higher  CD8+  immune  response  (1607  cells / /il-blood)  than  the  EQ2  when 
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ax  =  8  x  10  3  day  \  it  has  degraded  with  respect  to  the  CD4+  T-cell  count  (495.0  vs.  915.6 
cells//il-blood). 

If  we  now  start  simulations  with  ax  =  8  x  10-2  day-1  and  the  initial  conditions  (0,  0, 1400,  0, 
z/,  0,  0.01,  0),  then  the  solution  of  model  (2.1)  also  converges  to  either  uninfected  equilibrium 
or  “unhealthy”  infected  equilibrium  depending  on  v.  But  in  this  case,  the  solution  converges 
to  its  corresponding  “unhealthy”  infected  equilibrium  with  v  around  100.3  copies/ml-plasma 
(see  Fig.  2).  Hence,  the  viral  load  necessary  for  infection  depends  on  the  maximum  activation 
rate  ax',  low  values  of  «T  require  a  larger  viral  load  to  reach  an  infected  state.  Furthermore, 
when  infected,  smaller  values  of  ax  produce  “healthier”  infections  than  infections  with  larger 
ax  values,  in  terms  of  CD4+  T-cell  counts  and  viral  load  level. 

If  we  set  the  maximum  activation  rate  of  resting  CD4+  T-cells  lower,  such  as  ax  —  lx  10-3 
day-1,  and  keep  the  values  of  all  the  other  parameters  to  be  the  same  as  those  illustrated 
in  Table  2,  then  we  find  that  we  have  only  two  locally  a.s.  off-treatment  equilibria:  the 
uninfected  steady  state  (EQi)  and  “healthy”  infected  steady  state  (EQ2),  where  again  the 
uninfected  steady  state  is  the  same  as  with  ax  =  8  x  1CT3  day-1  (Table  3).  Even  though 
the  “healthy”  infected  steady  state  has  a  higher  viral  load  level  (310  copies/ml-plasma)  and 
lower  immune  response  (154  cells///l-blood)  than  it  does  when  ax  —  8  x  10-3  day-1,  the 
immune  response  still  controls  the  viral  load  to  maintain  it  below  a  detection  limit  of  400 
copies/ml-plasma.  This  equilibrium  also  has  much  higher  CD4+  T-cell  counts  (1290  cells / yul- 
blood)  than  the  case  when  ax  =  8  x  10-3  day-1.  If  we  start  the  simulation  with  the  initial 
conditions  (0,  0, 1400,  0,  z/,  0,  0.01,  0),  then  the  solution  of  model  (2.1)  converges  either  to  its 
uninfected  equilibrium  or  “healthy”  infected  equilibrium  based  on  the  value  of  z/,  and  it 
converges  to  this  “healthy”  infected  equilibrium  with  u  around  1.527  x  10'  copies/ml-plasma 
(Fig.  2). 

The  examples  above  demonstrate  the  existence  of  multiple  off-treatment  equilibria  and  il¬ 
lustrate  that  changing  the  value  of  a  parameter,  such  as  the  maximum  activation  rate  ax, 
has  an  effect  on  both  the  number  and  “health”  characterization  of  the  equilibria.  In  ad¬ 
dition,  the  initial  conditions  of  the  system  determine  which  equilibrium  is  attained  after 
initial  infection.  Across  a  population,  the  parameter  values  can  be  expected  to  vary  to  rep¬ 
resent  different  host  factors  and  host-virus  interaction  rates.  In  the  analysis  of  parameter 
ax,  we  End  that  a  person  with  a  lower  ax  value  requires  a  larger  viral  load  in  order  to  get 
infected  and,  once  infected,  attains  a  “healthier”  set  point  outcome,  in  terms  of  their  CD4+ 
T-cell  counts  and  viral  load  level,  than  those  with  higher  values  of  ax-  Analysis  with  other 
parameters,  such  as  the  expansion  factor  px,  reveals  similar  behavior. 

Multiple  equilibria  in  an  unhealthy  person  (a a  ^  0).  In  order  to  investigate  whether 
the  activation  due  to  a  non-HIV  antigen  can  affect  the  number  of  physical  equilibria  and 
their  local  asymptotic  stability,  we  performed  simulation  with  a  a  =  10-5  day-1  and  the 
values  of  all  the  other  parameters  as  specified  in  Table  2.  All  the  locally  a.s.,  physical,  off- 
treatment  equilibria  are  tabulated  in  Table  4,  which  indicates  that  we  still  have  three  locally 
a.s.  equilibria:  an  uninfected  steady  state  EQ i,  a  “healthy”  infected  steady  state  EQ-2  and 
an  “unhealthy”  infected  steady  state  EQ 3.  We  also  observe  that  EQ  1  has  a  nonzero  E2  even 
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though  bE2  <  $E2  and  a  nonzero  T\,  which  is  not  observed  for  equilibrium  EQi  when  oa  =  0 
(Table  3).  Otherwise,  the  equilibria  for  the  two  cases  (a^  =  8  x  10-3  day-1  and  a  a  =  0  or 
1  x  10-5  day-1)  are  very  similar. 


state 

aA 

1  x  10-5 

1  x  10-2 

EQ i 

eq2 

eq3 

EQ i 

eq2 

EQ3 

T\ 

0.6710 

184.5 

13.13 

266.7 

277.2 

13.12 

T* 

0 

0.05621 

5.332 

0 

0.05592 

6.567 

t2 

1398 

730.6 

424.0 

555.6 

404.2 

228.2 

t; 

0 

0.04388 

2.982 

0 

0.01835 

1.777 

Vj 

0 

265.0 

28420 

0 

248.2 

35010 

Vni 

0 

0 

0 

0 

0 

0 

E1 

0.009938 

799.9 

0.04573 

0.009121 

1160 

0.04113 

e2 

0.001562 

98.41 

0.002862 

0.02198 

147.6 

0.002630 

local  stability 

a.s. 

a.s. 

a.s. 

unstable 

a.s. 

a.s. 

Table  4:  Off-treatment  steady  states  for  model  (2.1)  with  =  1  x  10-5  and  1  x  10-2  day-1 
and  the  values  of  all  the  other  parameters  as  specified  in  Table  2.  Non-physical  and  unstable 
steady  states  are  omitted,  except  in  the  case  of  EQi  when  cza  =  1  x  1CT2  day-1,  which  is 
unstable.  Note  that  ar  =  8  x  10-3  day-1  for  this  analysis. 

If  we  start  simulations  with  the  initial  conditions  (0.6710,  0, 1398,  0,  is,  0, 0.009938,  0.001562), 
then  the  solution  of  model  (2.1)  converges  either  to  its  uninfected  equilibrium  EQi  or  “un¬ 
healthy”  infected  equilibrium  EQ 3  based  on  the  value  of  is,  and  it  will  not  converge  to  EQ 3 
until  is  is  close  to  4000  (Fig.  3),  as  compared  to  6576  copies/ml-plasma  when  a  a  =  0.  Hence, 
a  person  infected  with  a  non-HIV  virus  becomes  HIV-infected  at  lower  viral  loads  than  when 
no  other  infection  is  present. 

If  we  take  the  value  of  a  a  larger,  such  as  a  a  =  1  x  10-2  day-1,  then  the  uninfected  equilibrium 
is  no  longer  stable,  but  we  still  have  the  locally  a.s.  “healthy”  and  “unhealthy”  steady  states 
designated  as  EQ2  and  EQ3  in  Table  4,  respectively.  If  we  start  simulations  with  the 
initial  conditions  (266.7,  0, 555.6, 0,  is,  0,  0.009121, 0.02198),  then  the  solution  of  model  (2.1) 
converges  to  EQ3,  regardless  of  the  initial  viral  load  is  (Fig.  3).  This  means  that  when  the 
system  is  in  the  state  EQ i,  introduction  of  even  the  smallest  amount  of  virus  will  cause  the 
system  to  converge  to  the  “unhealthy”  equilibrium  EQ3. 

Furthermore,  Table  4  indicates  that  as  the  activation  rate  by  non-HIV  antigen  a  a  becomes 
larger,  the  locally  a.s.  “healthy”  steady  state  has  degraded,  in  terms  of  the  total  CD4+ 
T-cells  counts  (915.1  vs.  681.1  cells//d-blood),  while  the  viral  load  has  improved  (265.0  vs. 
248.2  copies/ml-plasma),  as  well  as  the  immune  response  (898.4  vs.  1307  cells//il-blood).  In 
addition,  the  locally  stable  “unhealthy”  steady  state  has  worsened,  in  terms  of  a  much  lower 
CD4+  T-cell  count  (445.4  vs.  249.4  cells/pl-blood),  much  higher  viral  load  set  point  (28420 
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Figure  3:  Phase  diagram  showing  equilibrium  attained  as  a  function  of  the  initial  viral  load 
v  and  the  parameter  am  Initial  conditions  for  simulations,  corresponding  to  the  uninfected 
equilibrium  EQi,  are  calculated  for  each  value  of  am  Note  that  clt  =  8  x  10-3  day-1  for 
these  simulations. 

vs.  35010  copies/ml-blood),  and  much  lower  immune  response  (0.04859  vs.  0.04376  cells / yul- 
blood).  Therefore,  at  this  point  increasing  a  a  leads  to  similar  behavior  as  does  increasing 
ax  when  a  a  =  0. 

The  examples  above  summarize  an  investigation  as  to  the  effect  of  a  non-HIV  infection 
(. cla  0)  on  the  outcome  of  an  HIV-infection.  We  find  that  an  unhealthy  person  (a^  ^  0) 
requires  a  smaller  initial  viral  load  v  to  sustain  an  HIV  infection  and,  once  infected,  attains  a 
worse  outcome,  in  terms  of  CD4+  T-cell  counts  and  viral  load  levels,  than  a  healthy  person 
{a a  =  0).  We  find  that  multiple  equilibria,  including  an  uninfected  equilibrium,  exist  in  the 
presence  of  a  small  non-HIV  infection  (a^  =  1  x  1CT5  day-1).  However,  for  larger  levels  of 
non-HIV  infection  (a a  =  1  x  10-2  day-1),  the  uninfected  equilibrium  is  unstable  and  even 
the  smallest  initial  viral  load  will  lead  to  an  “unhealthy”  equilibrium  state  (Fig.  3). 


4  Viral  Blips 


Adherence  to  a  regimen  of  HAART  suppresses  the  viral  loads  of  most  infected  HIV  patients 
below  the  level  of  detection  (<400  or  <50  copies/ml-plasma  depending  on  the  assay  used) 
by  standard  assays  .  However,  a  number  of  these  well-suppressed  patients  experience  unex¬ 
plained  “viral  blips”  while  on  therapy  [15,  20].  In  a  study  [15]  of  123  patients,  these  viral 
blips  are  estimated  to  have  a  duration  of  approximately  two  to  three  weeks  with  a  mean  blip 
amplitude  of  158  ±132  copies/ml-plasma.  Furthermore,  it  is  observed  that  the  blip  frequency 
inversely  correlates  with  CD4+  T-ccll  counts.  In  this  section,  we  investigate  how  infection 
with  a  non-HIV  antigen  can  lead  to  viral  blips  for  those  patients  who  are  on  treatment  and 
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have  successfully  suppressed  their  viral  loads  to  undetectable  levels.  To  do  so,  we  use  the 
states  of  a  locally  a.s.  “healthy”  infected  equilibrium  (201,  0.056,  730,  0.040,  288,  0,  240,  29.2), 
corresponding  to  =  0.7,  e2  =  0,  =  0,  and  the  values  of  all  the  other  parameters  as 

specified  in  Table  2,  as  the  initial  conditions  for  simulations  in  which  we  allow  parameter 
cla  to  be  a  function  of  time  t,  while  keeping  all  other  parameter  values  to  be  the  same.  As 
we  drive  the  system  (2.1)  with  the  non-HIV  infection  cia(£)  we  monitor  the  states  of  the 
system,  particularly  the  viral  load.  While  the  cause  of  viral  blips  has  not  yet  been  resolved, 
one  proposed  mechanism  [20]  posits  that  viral  blips  could  be  due  to  an  increase  in  activated 
CD4+  cells  as  a  result  of  secondary  infection  or  vaccination.  Our  approach,  introducing 
a  time-dependent  non-HIV  infection  ayi(f)  is  consistent  with  this  proposed  mechanism  and 
demonstrates  that  our  model  supports  such  a  scenario. 

Viral  blip  caused  by  a  single  non-HIV  infection.  Figure  4(a)  illustrates  the  case 
of  a  one-time  non-HIV  infection  occurring  on  days  20  through  50,  with  a  peak  value  of 
0.003  day'1  (days  30  to  40).  Figure  4(b)  depicts  the  HIV  viral  progression  Vj(t)  resulting 
from  this  infection.  In  Fig.  4(b)  we  see  that  there  is  a  small  delay  before  the  viral  load 
begins  to  increase.  This  delay  is  reasonable  since  there  is  no  activation  term  (a^)  in  the  Vj 
compartment  of  model  (2.1)  and  the  viral  load  can  only  increase  after  T*  increases.  About 
10  days  after  the  start  of  the  infection  (day  30),  the  viral  load  rises  above  the  detection  limit 
(400  copies/ml-plasma),  reaching  its  peak  value  15  days  post-infection  (day  35).  This  viral 
blip  drops  below  the  detection  limit  20  days  post-infection  (day  40),  before  the  infection  has 
completely  cleared. 


Figure  4:  (a)  Activation  rate  a  Ait)  by  a  non-HIV  antigen.  The  duration  of  the  infection  is 
30  days;  (b)  Viral  load  Vj{t)  (solid  line),  censored  data  level  (horizontal  dashed  line),  and 
infection  start  and  stop  times  (vertical  dash-dot  lines). 


To  further  investigate  the  dynamics  of  the  viral  blip,  we  plot  the  other  model  compartments 
for  the  example  shown  in  Fig.  4  (we  omit  compartment  Vw,  since,  with  e2  =  0  and  Vv/(0)  = 
0,  Vnj  =  0).  In  Fig.  5,  it  can  be  seen  that  the  T-ccll  compartments  (T),  7\*,  T2,  and  T2*) 
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respond  at  the  start  of  the  infection  (day  20).  This  is  not  surprising  because  the  non-HIV 
antigen  activation  term  a  a  appears  in  the  dynamical  equations  for  these  compartments  in 
model  (2.1).  It  is  also  reasonable  that  the  resting  (naive  and  memory)  T-cell  compartments 
(T2  and  T2* )  initially  decrease,  while  the  activated  T-cell  compartments  (7\  and  Tj*)  increase, 
since,  at  the  start  of  the  infection,  resting  cells  are  reactivated  (memory)  to  fight  the  infection 
or  activated  through  infection  (naive).  There  is  also  a  short  delay  between  the  peak  value 
of  T2  and  the  peak  value  of  Tj*  which  is  likely  due  to  the  time  delay  in  the  differentiation  of 
Tj*  into  T2* .  Comparing  compartment  Tj*  of  Fig.  5  with  compartment  Vj  in  Fig.  4,  we  find 
that  these  two  compartments  exhibit  very  similar  time-dependent  behavior.  Again,  this  is 
to  be  expected  because  the  only  compartment  that  produces  virus  is  Tj*. 

If  we  focus  our  attention  on  the  uninfected  T-cell  compartments  (T\  and  T2  in  Fig.  5), 
we  see  that  the  peak  of  T\  corresponds  to  the  valley  of  T2.  The  fact  that  Tj  increases 
and  T2  decreases  for  most  of  the  infection  implies  that  the  effects  of  the  activation  terms 
Pt  (otVi/(Vi  +  Ky)  +  a  a)  T2  and  (cltVi/ (Vj  +  Kv)  +  a  a)  T2  dominate  over  the  differentia¬ 
tion  terms  JtTi  and  other  loss  terms  in  model  (2.1)  during  the  infection.  That  is,  the  T2 
compartment  is  losing  more  through  activation  then  it  is  gaining  through  differentiation  of 
Tj,  and  vice  versa  for  T\.  The  situation  is  reversed  as  the  infection  clears  (day  50),  where 
there  is  a  slow  decay  of  T\  and  corresponding  rise  in  T-2  back  toward  the  equilibrium  values, 
as  the  decay  of  Ti  contributes  to  the  rise  of  T2  through  differentiation. 
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Figure  5:  Model  dynamics  with  t\  =  0.7,  e2  =  0,  and  the  values  of  all  the  other  parameters 
as  specified  in  Table  2.  Activation  rate  by  a  non-HIV  antigen  a^(t)  is  as  depicted  in  Fig. 
4(a).  Vertical  dashed  lines  indicate  the  start  and  stop  times  of  the  non-HIV  infection. 
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It  is  interesting  to  contrast  the  above  behavior  of  T\  and  T2  with  the  behavior  of  their  infected 
counterparts  (Tj"  and  T2).  Unlike  the  T2  compartment,  which  decreases  throughout  most  of 
the  infection,  the  latently-infected  T2  compartment  only  decreases  for  a  brief  period,  then 
rises  and  peaks  near  the  midpoint  of  the  infection  (Fig.  5).  In  behavior  similar  to  that  in 
the  T2  compartment,  the  initial  drop  in  T2  is  due  to  activation  of  naive  cells  and  reactivation 
of  memory  cells,  although,  in  fact,  the  T2  compartment  mostly  represents  infected  memory 
cells,  since  the  rate  of  infection  of  naive  cells  is  very  low.  This  activation  is  a  source  term 
for  the  Tj"  compartment  and  a  loss  term  for  the  T2*  compartment,  but  there  is  also  an 
important  source  term  for  T2*  that  begins  to  dominate  over  the  loss  due  to  activation.  As 
the  T\  compartment  increases,  there  is  an  increased  number  of  T\  cells  becoming  infected 
and  adding  to  the  Tj"  compartment,  which,  in  turn,  leads  to  increased  differentiation  to  the 
T2*  compartment.  The  net  result  is,  that  after  a  short  decrease,  the  T2  cells  increase,  even  as 
Tj"  continues  to  increase,  as  the  differentiation  (source)  term  7 ^Tj"  in  the  T2  compartment 
dominates  over  the  activation  (loss)  term  (o,tVi/(Vi  +  Kv)  +  a  a)  T2  in  model  (2.1).  The  net 
result  is  that  the  latently-infected  T2  compartment  is  supplemented,  not  depleted  during  the 
secondary  infection.  However,  as  a  a  drops,  the  loss  terms  begin  to  dominate  in  the  Tj"  and 
T2  compartments  and  both  drop  below  their  equilibrium  values  before  slowly  rising  back  to 
the  equilibrium  values.  Indeed,  150  days  after  the  end  of  the  infection  only  Tj"  has  returned 
to  its  equilibrium  value  (Fig.’s  4  and  5). 

We  now  focus  on  the  CD8+  compartments  (El  and  E2),  which  are  responsible  for  suppressing 
the  viral  blip.  In  Fig.  5  we  can  see  that  E\  and  E2  respond  to  the  infection  with  longer  delays 
than  the  other  compartments.  This  is  reasonable  since  the  activation  of  the  CD8+  memory 
cells  ( E2 )  does  not  depend  on  the  non-HIV  antigen  cla  (these  are  HIV-specific  memory  cells) 
and  activation  will  occur  in  response  to  a  change  in  Vj.  The  CD8+  compartments  also  have 
a  more  complicated  relationship  to  T\  and  Tj".  In  Fig.  5,  we  see  that  both  compartments 
increase  through  most  of  the  infection,  although  E2  has  a  short  period  in  which  it  decreases 
at  the  beginning  of  the  infection.  This  decrease  is  probably  due  to  the  loss  of  cells  due  to 
the  activation  term  aEViE2/{Vi  +  Ky)  in  model  (2.1).  The  subsequent  increase  in  both 
CD8+  compartments  indicates  that  the  additional  source  term  bpAT'lEl  j (Tj"  +  Kbl)  in  the 
Ei  compartment  and  the  resulting  differentiation  (rjEEi(Ti  +  Tj")/(Tx  +  Tj"  +  A7))  to  E2 
dominate  over  the  activation  (loss)  term  aEVjE2/ {Vj  +  Kv)  in  the  E2  compartment.  It  is 
interesting  to  note  that  the  peak  value  of  E-2  occurs  after  the  infection  has  cleared  (around 
day  60).  As  the  viral  load  drops,  the  loss  in  E2  due  to  activation  decreases,  while  the  higher 
levels  of  Ei  continue  to  provide  a  source  for  E2  through  differentiation. 

The  viral  blips  produced  with  this  model  last  about  10  days,  not  2-3  weeks  as  estimated  in 
the  literature  [15].  Some  of  this  discrepancy  may  be  due  to  the  fact  that  we  are  using  a  higher 
detection  limit  (400  RNA  copies/ml-plasma)  than  the  experimental  studies.  Also,  as  we  have 
seen,  the  dynamics  of  the  V}  compartment  closely  follow  those  of  the  Tj"  compartment,  which 
is  the  only  source  for  Vj.  This  may  explain  why  the  the  sharply  rising  viral  load  in  Fig.  4 
is  not  followed  by  a  slower,  two-phase  decay  as  observed  in  [20],  where  there  is  another 
compartment  of  less-rapidly  changing  “chronic”  cells  that  also  produce  virus.  Also,  in  an 
effort  to  facilitate  fitting  to  clinical  data,  we  have  combined  naive,  memory,  activated,  and 
helper  (uninfected)  T-cells  in  two  compartments  (Ti,  T2).  Perhaps  further  subdividing  the 
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CD4+  T-ccll  compartments  or  adding  more  compartments,  such  as  a  chronically  infected 
T-cells  pool  or  antigen-specific  cells  and  non-antigen  specific  cells  [20],  would  remedy  this 
discrepancy. 

Viral  blips  caused  by  sequential  non-HIV  infections.  We  next  investigate  the  effect 
of  sequential  non-HIV  infections  by  simulating  two  such  back-to-back  infections.  The  first 
infection  occurs  from  days  20  through  50  and  stays  at  its  peak  value  (0.003  day-1)  from  days 
30  through  40,  while  the  second  infection  occurs  from  days  80  through  110  and  stays  at  its 
peak  value  (0.003  day-1)  from  days  90  through  100  (Fig.  6(a)).  As  seen  in  Fig.  6(b),  the 
first  infection  leads  to  a  viral  blip  at  exactly  the  same  time  as  that  in  Fig.  4  and  the  second 
infection  leads  to  a  viral  blip  (starting  at  day  95)  with  the  same  duration  as  the  first,  but 
with  a  smaller  amplitude. 


(a)  (b) 


Figure  6:  (a)  Activation  rate  by  non-HIV  antigens  a^(f).  The  duration  of  each  infection  is 
30  days,  with  the  first  beginning  on  day  20  and  the  second  on  day  80;  (b)  Viral  load  Vj{t) 
(solid  line)  and  censor  data  level  (horizontal  dashed  line). 


To  investigate  whether  the  short  time  interval  between  these  two  infections  leads  to  the 
suppressed  amplitude  of  the  second  viral  blip,  we  also  considered  a  case  where  the  first 
infection  is  the  same  as  that  in  Fig.  6  (days  20  through  50)  but  the  second  infection  occurs 
much  later  (days  200  through  230).  The  simulation  results  are  illustrated  in  Fig.  7,  which 
indicates  that  the  amplitude  of  the  second  viral  blip  now  is  slightly  larger  than  the  first  viral 
blip. 
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(a)  (b) 


Figure  7:  (a)  Activation  rate  by  non-HIV  antigens  a  Ait).  The  duration  of  each  infection  is 
30  days,  with  the  first  beginning  on  day  20  and  the  second  on  day  200;  (b)  Viral  load  Vi(t) 
(solid  line)  and  censor  data  level  (horizontal  dashed  line). 


Based  on  these  examples,  we  see  that  the  frequency  of  non-HIV  infections  affects  the  fre¬ 
quency  of  the  viral  blips  and  that  the  time  interval  between  the  two  infections  can  affect 
the  amplitude  of  the  second  viral  blip.  In  particular,  the  shorter  the  time  interval  between 
two  infections,  the  smaller  the  amplitude  of  the  second  viral  blip.  This  suggests  that  it  is 
possible  that  there  may  be  only  one  viral  blip  detected  if  the  time  interval  between  the  two 
infections  is  short  enough. 

To  further  investigate  the  dynamics  of  two  sequential  infections,  we  plot  the  other  model 
compartments  for  the  example  shown  in  Fig.  7.  In  Fig.  8  we  observe  that  the  E\  and 
E2  compartments  are  still  below  their  equilibrium  values  when  the  second  infection  occurs. 
This  may  explain  the  increased  amplitude  of  the  second  viral  blip.  In  the  case  of  the  more 
closely-spaced  second  infection  of  Fig.  6(a),  the  start  of  which  is  indicated  with  the  dotted 
vertical  line  in  Fig.  8,  we  can  see  that  E\  and  E2  are  above  their  equilibrium  levels.  This 
may  explain  the  decreased  amplitude  of  the  second  viral  blip  in  Fig.  6(b). 
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Figure  8:  Model  dynamics  with  t\  =  0.7,  62  =  0,  and  the  values  of  all  the  other  parameters  as 
specified  in  Table  2.  Activation  rate  by  a  non-HIV  antigen  a^f)  is  as  depicted  in  Fig.  7(a). 
Vertical  dashed  lines  indicate  the  start  and  stop  times  of  the  non-HIV  infections  for  this 
example.  The  vertical  dotted  line  indicates  the  start  time  (80  days)  of  the  second  infection 
in  Fig.  6. 
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The  effect  of  the  form  of  the  non-HIV  infection  on  the  viral  blip.  As  one  might 
expect,  both  the  shape  and  the  peak  value  of  a  Ait)  affect  the  amplitude  and  duration  of  the 
viral  blip.  As  depicted  in  Fig.  9,  we  find  that  with  the  same  shape  a^(t)  but  with  different 
amplitude  (a^  and  a®),  the  viral  blip  tends  to  have  higher  amplitude  and  a  little  bit  wider 
duration  as  peak  value  of  a^(t)  gets  higher  (see  Vj1'1  and  vf:i>  in  Fig.  9(b)).  Furthermore, 
we  see  that  amplitude  of  the  viral  blip  is  lower  as  the  non-HIV  infection  becomes  longer  and 
the  peak  value  of  the  viral  blip  also  occurs  earlier  (vjl>  and  Vj2'1  in  Fig.  9(b)). 


Figure  9:  (a)  Activation  rate  by  non-HIV  antigen  <za(£);  (b)Viral  load  Vi(t). 


CD4+  help  during  a  viral  blip.  In  model  (2.1),  the  differentiation  of  Ei  cells  to  E2  mem- 

1 1  1  _|_/T>* 

ory  cells  occurs  through  the  term  1et  +T*+K  E\ ,  which  is  dependent  on  the  help  of  activated 
CD4+  cells  (T\  and  T*).  In  the  case  that  T\  and  T*  are  zero,  there  is  no  differentiation  of 
the  effector  Ei  cells  into  memory  E2  cells.  On  the  other  hand,  in  the  limit  of  large  T\  +  T£, 
the  maximum  rate  of  differentiation  is  attained.  In  this  section  we  explore  what  effect 
the  CD4+  help  has  on  the  CD8+  immune  response  and  we  do  so  in  the  context  of  a  viral 
blip  induced  by  a  non-HIV  infection  a  Ait)-  We  use  the  parameter  K1  to  modify  the  effects  of 
the  CD4+  help.  The  case  of  JF7  =  0  simulates  a  condition  where  CD4+  help  is  not  required 
for  differentiation  to  CD8+  memory,  while  the  case  of  K1  very  large  simulates  a  condition 
where  differentiation  to  CD8+  memory  is  impaired  due  to  lack  of  CD4+  help.  Thus,  we  use 
the  states  of  locally  a.s.  “healthy”  infected  equilibria,  corresponding  to  =  0.7,  e2  =  0, 
cla  =  0,  various  values  of  Jl7,  and  the  values  of  all  other  parameters  as  specified  in  Table  2, 
as  the  initial  conditions  for  simulations  in  which  we  allow  the  parameter  a  a  to  be  a  function 
of  time  t.  As  we  drive  the  system  (2.1)  with  the  non-HIV  infection  a  Ait)  we  monitor  the 
states  of  the  system. 

In  Fig.  10(a)  we  plot  the  “healthy”  infected  (EQ2)  equilibrium  values  of  the  Vj,  Ei,  and 
E2  compartments  as  a  function  of  K1.  As  expected,  as  K1  increases,  corresponding  to 
increasing  impairment  of  CD8+  differentiation  to  memory,  the  equilibrium  value  of  the  E2 
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compartment  decreases.  However,  the  overall  effect  on  E\  is  negligible  with  only  a  slight 
rise  as  K1  increases.  Despite  the  fact  that  Ei  stays  relatively  constant,  we  see  that,  as  K1 
increases,  the  equilibrium  value  of  the  viral  load  increases  until  reaching  a  plateau  value  at 
K-,  «  1  x  105  cells//il-blood.  Thus,  even  though  the  number  of  effector  cells,  which  directly 
participate  in  the  removal  of  infected  T-cells,  does  not  decrease  as  the  CD8+  differentiation 
to  memory  is  impaired,  the  overall  viral  load  increases. 


(a)  (b) 

Figure  10:  (a)  The  effect  of  Ky  on  the  “healthy”  infected  (EQ2)  values  of  the  Vj  (circles), 
Ei  (crosses),  and  E-2  (asterisks)  compartments.  Equilibrium  values  correspond  to  the  case 
c±  =  0.7,  e2  =  0,  a  a  =  0,  and  the  values  of  all  other  parameters  as  specified  in  Table  5;  (b) 
Relative  change  of  Vj  (circles),  E\  (crosses),  and  E-2  (asterisks)  compartments  as  a  function 
of  a;  . 


In  Fig.  10(b)  we  plot  the  relative  changes  of  the  various  compartments  as  a  function  of  Ji7. 
Relative  change  is  defined  as  (peak  value  minus  initial  value)/ (initial  value).  In  this  figure 
we  can  see  that,  as  K1  increases  the  relative  change  of  the  peak  viral  load  increases.  The 
relative  change  of  the  E2  memory  compartment  also  increases  (Fig.  10(b)),  but  the  infected 
equilibrium  value  was  quite  small  to  begin  with  (Fig.  10(a)).  The  relative  change  in  the 
effector  Ei  compartment  decreases  as  Ji7  increases.  The  results  of  this  analysis  demonstrates 
that  in  our  model,  impaired  differentiation  to  CD8+  memory  results  in  a  degraded  “healthy” 
equilibrium  and  leads  to  a  larger  viral  blip  following  a  secondary  infection. 
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5  Fitting  the  Model  to  Clinical  Data 


When  proposing  a  new  model,  a  necessary  but  often  difficult  requirement  is  validation  of 
the  model  with  clinical  data.  One  way  to  do  this  is  to  use  the  first  part  of  the  available 
longitudinal  data  for  a  given  patient  to  estimate  patient  specific  parameters.  One  then  uses 
these  parameters  in  model  simulations  to  accurately  predict  the  subsequent  progression  of 
the  disease  as  represented  in  the  remainder  of  the  longitudinal  data  for  that  patient.  In  this 
section  we  test  the  capability  of  model  (2.1)  to  fit  the  clinical  data  of  patients  and  to  predict 
the  viral  progression.  Of  course,  the  dynamics  of  HIV  progression  may  totally  change  as  time 
evolves  (  e.g.,  such  as  transfer  from  “unhealthy”  steady  state  region  of  attraction  to  one  for 
a  “healthy”  state).  In  this  section,  we  test  the  predictive  capabilities  of  the  model  with 
parameters  that  have  been  estimated  from  only  half  of  the  longitudinal  data  by  comparing 
simulations  to  the  full  set  of  clinical  observations. 

To  obtain  patient  specific-parameter  estimates  for  the  model,  we  used  individual  patient 
data  (either  partial  or  full  longitudinal  data)  and  carried  out  an  inverse  problem.  The 
data  for  our  investigations  come  from  Massachusetts  General  Hospital  (MGH),  where  all 
the  patients  enrolled  in  the  study  are  symptomatic  with  acute  or  early  HIV-infection  (for 
more  detailed  information  of  these  data,  the  interested  reader  is  referred  to  [3,  21,  27]).  In 
summary,  nearly  all  subjects  in  the  study  undergo  combination  therapy  and  many  have  at 
least  one  treatment  interruption.  The  available  clinical  data  include  total  CD4+  T-cell  count 
and  total  RNA  copies,  where  for  model  (2.1)  the  total  CD4+  T-ccll  counts  are  represented 
by  zi(t]q)  =  Ti(f;g)  +  Tj*(f;g)  +  T2{t\q )  +  Tf(t\q)  and  total  RNA  copies  are  represented 
by  z2(t:  q)  =  V}(f;g)  +  VNi(t',q).  If  the  measurements  of  RNA  copies  are  below  the  limit 
of  quantification  for  the  assay  used  (400  copies/ml-plasma  for  a  standard  assay  and  50 
copies/ml-plasma  for  an  ultra-sensitive  assay),  then  the  observed  viral  load  value  is  censored 
to  be  at  its  detection  limit;  that  is,  in  these  cases  the  observed  values  do  not  represent  the  true 
data  values  anymore.  Furthermore,  observations  of  viral  load  and  CD4+  may  not  be  at  the 
same  time  points  and  the  observation  times  and  intervals  vary  substantially  among  patients. 
So,  in  general,  for  patient  number  j  we  have  CD4+  T-cell  data  pairs  (tf ,  yf),  i  —  1,  •  •  •  ,  N(, 
and  potentially  different  time  point  viral  RNA  data  pairs  (fif  ,y2),  i  —  1,  •  •  •  ,  At).  Hence,  the 
clinical  data  for  carrying  out  the  inverse  problem  involves  partial  observations ,  measurements 
from  combined  compartments ,  and  highly  censored  viral  load  measurements. 

As  one  might  expect,  if  a  patient  does  not  have  a  sufficient  number  of  observations  or  does 
not  undergo  a  therapy  interruption  during  the  observation  period  used  to  fit  the  model,  then 
it  is  difficult  to  obtain  dynamically  dependent  parameters  to  use  to  predict  a  later  disease 
progression  involving  both  therapy  and  treatment  interruption.  Hence,  for  our  purposes 
we  focused  on  14  patients  wherein  each  patient  has  ten  or  more  each  of  CD4+  and  viral 
load  measurements,  and  have  at  least  one  on/off  treatment  schedule  in  the  first  half  of  their 
longitudinal  data.  Moreover,  for  each  of  these  patients,  N(  is  not  substantially  different 
than  N%.  The  technique  we  use  in  this  paper  is  adapted  from  the  one  in  [3],  where  the 
authors  developed  a  statistically-based  censored  data  method  (an  expectation  maximization 
algorithm)  combined  with  an  ordinary  nonlinear  least-squares  technique.  We  note  that 
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the  variance  a\  in  CD4  measurements  and  the  variance  a2  in  viral  load  measurements  are 
likely  to  be  different  due  to  assay  differences.  Hence,  in  this  paper  we  use  the  expectation 
maximization  algorithm  based  on  Maximum  Likelihood  Estimation  for  (q,af,a2)  which, 
under  normality  assumptions  on  the  errors,  results  in  a  weighted  least-squares  technique 
(see  [30])  with  solution  given  by 

n{  n32 

A  E  «)l2  +  4  E  l»?  -  *(*?;?)  I2 

07  2 — '  07.  2 — ' 

1  i=l  2  i=l 

iV£ 

i=l 

for  the  logio-transformed  system  of  model  (2.1)  for  patient  j,  where  y{  =  log 10yjJ  and 
zi =  logio(^i i  =  I,’  ’*'  )A^,  and  =  log10$f|  z2(f2J;q)  =  log10(z2(t^ ;  g)), 
i  =  1,  •  •  •  ,  jVg.  As  noted  in  [3],  by  using  a  log-transformed  system  one  can  resolve  a  problem 
of  states  becoming  unrealistically  negative  due  to  round-off  error:  nonnegative  solutions 
of  this  model  should  stay  so  throughout  numerical  simulation.  This  approach  also  enables 
efficient  handling  of  unrealistic  cases  where  states  get  infinitesimally  small  during  integration 
due  to  parameters  selected  by  optimization  algorithms.  From  a  statistical  point  of  view,  log 
transformation  is  a  standard  technique  to  render  the  observations  more  nearly  normally 
distributed,  which  also  supports  use  of  the  weighted  least  squares  criterion  as  an  equivalent 
to  maximum  likelihood  estimation. 

The  expectation  maximization  (EM)  algorithm  is  outlined  below.  To  simplify  the  notation, 
we  drop  the  patient  index  j  in  this  algorithm  description.  The  following  notation  will  be  used 
in  the  algorithm:  the  relevant  censoring  point  at  time  f  is  represented  by  Ll  and  y*  is  the 
indicator  function  for  the  set  {y2  >  Lr},  <fi  denotes  the  standard  normal  probability  density 
function  and  <f>  is  the  corresponding  cumulative  distribution  function.  For  each  patient,  we 
carried  out  the  following  parameter  estimation  algorithm: 


(5.1) 


qJ  =  argmin 
q£Q 


•  (Step  1)  Create  adjusted  data  yl  by  replacing  censored  y2  values  by  V / 2,  and  use 
ordinary  least  squares  to  estimate  q ^  using  both  CD4+  data  y\  and  adjusted  viral 
RNA  data  yl  (which  now  includes  replaced  censored  values). 

Ni  N2 
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Obtain  an  initial  estimate  for  of  and  o\  from 
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respectively.  Set  k  =  0. 
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•  (Step  2)  Define  =  Z2(tl2',q^)  and  =  L  „(% — •  Update  the  data  and  residuals 
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•  (Step  3)  Update  the  estimate  of  g  to  g(fc+1^  by  performing  the  weighted  least  squares 
minimization  in  the  parameters  q 


g(fc+1)  =  argmin 
q&Q 
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and  computing 
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If  relative  changes  in  q,  &i  and  are  small,  terminate.  Otherwise  set  k  —  k  +  1  and 
then  go  to  Step  2. 


For  more  details  about  the  EM  algorithm  and  how  to  carry  it  out,  interested  readers  are 
referred  to  [3,  13,  23]. 

Note  that  model  (2.1)  has  31  model  parameters  and  8  initial  conditions.  The  total  number 
of  measurements  (sum  of  number  of  measurements  of  CD4+  and  number  of  measurements 
of  RNA  copies)  in  the  first  half  of  the  longitudinal  data  varies  from  42  to  154  in  these  14 
patients,  which  means  it  is  difficult  to  estimate  all  39  of  these  parameter  values  for  some  of 
these  patients  by  using  half  of  the  longitudinal  data  set.  Hence,  we  first  try  to  estimate  all  the 
31  model  parameters  and  8  initial  conditions  for  each  of  the  14  patients  by  applying  the  EM 
algorithm  to  the  full  longitudinal  data  set.  We  then  try  to  fix  some  model  parameters  and 
initial  conditions  at  the  population  averages  across  these  patients,  and  attempt  to  estimate 
all  the  remaining  model  parameters  and  initial  conditions  for  each  patient  by  applying  the 
EM  algorithm  separately  to  each  of  the  first  half  and  the  full  longitudinal  data  set.  Note 
that  there  exists  biological  variation  in  all  parameters  across  the  patients,  and  there  also 
exist  high  correlations  among  some  of  these  parameters  such  as  the  RNA  copies  produced 
per  infected  cells  Nt  and  the  virus  natural  death  rate  c.  We  also  observe  that  sensitivity 
with  respect  to  some  of  these  parameters  may  be  highly  time-dependent.  For  example,  the 
dynamics  of  the  model  is  much  more  sensitive  to  the  treatment  efficacies  and  62  in  the 
treatment  periods  than  it  is  in  the  off-treatment  periods.  All  of  these  considerations  make 
it  difficult  to  choose  a  priori  which  parameters  can  be  fixed.  In  this  paper,  we  empirically 
chose  some  parameters  (such  as  the  saturation  parameters)  to  which  the  model  appears  to 
be  relatively  insensitive  to  take  as  fixed.  Table  5  specifies  all  the  fixed  parameters  and  their 
corresponding  values.  After  we  obtain  the  two  sets  of  model  parameters  (corresponding  to 
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Tabic  5:  Average  model  parameter  values  (19)  and  initial  conditions  (4)  used  in  model  fitting 
with  half  and  full  longitudinal  data  set. 

use  of  the  first  half  of  the  data  and  the  full  set  of  data,  respectively,  in  inverse  problems), 
we  simulate  the  trajectory  over  the  full  time  span  of  the  patient’s  observations  by  using  the 
parameter  values  obtained  with  these  two  data  sets  and  compare  their  ability  to  describe 
the  experimental  results. 

Figures  11,  12,  and  Fig.’s  13-24  in  the  Appendix  illustrate  model  fits  for  all  14  patients.  These 
figures  reveal  that  both  the  fits  to  CD4+  T-cells  and  the  fits  to  viral  load  are  reasonable 
when  using  either  half  data  or  full  data  to  estimate  the  parameter  values.  However,  there  is  a 
small  amount  of  under-prediction  or  over-prediction  for  some  patients  when  using  parameters 
estimated  from  the  half  data  set.  For  example,  model  fits  obtained  by  using  half  and  full  time 
series  data  for  patient  2  are  shown  in  Fig.  11,  where  it  can  be  seen  that  by  using  parameters 
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Figure  11:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 


estimated  from  the  first  half  data  set,  the  predicted  values  for  viral  load  data  are  slightly 
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higher  and  the  predicated  values  for  CD4+  T-cell  are  slightly  lower  than  the  observed  data 
values  and  the  fitted  values  from  using  parameters  estimated  from  the  full  data,  even  though 
there  are  two  off-treatment  periods  in  the  first  half  of  the  data  set.  This  discrepancy  may  be 
due  to  the  fact  that  both  off-treatment  periods  are  short  and,  therefore,  provide  very  little 
information  about  the  off-treatment  behavior.  Note  that  the  values  of  most  observations  of 
viral  loads  in  the  off-treatment  periods  during  the  first  half  data  set  are  higher  than  those 
observed  in  the  second  half  data  set  and  this  makes  it  difficult  to  correctly  predict  future 
off-treatment  trends.  Hence,  it  is  not  just  the  number  of  off-treatment  periods  in  the  first 
half  data  set  but  also  the  length  of  each  off-treatment  period  and  the  number  of  observations 
in  these  periods  that  determines  the  accuracy  of  the  predictions  of  off-treatment  viral  load 
levels. 

Figure  12  compares  the  model  fit  and  prediction  obtained  by  using  the  half  data  set  to  the  fit 
obtained  from  parameter  estimation  with  full  data  set  for  patient  10.  Even  though  there  is 
only  one  off-treatment  period  in  the  first  half  data  set,  the  period  is  sufficiently  long  and  the 
data  set  sufficiently  rich  that  the  full  data  fit  and  half  data  predictions  of  both  CD4+  T-cell 
and  viral  load  agree  quite  well  in  the  second  half  data  set.  Hence,  the  results  demonstrated 
in  Fig.’s  11  and  12  suggest  that  in  order  to  have  accurate  predictive  capability  we  need  to 
obtain  a  sufficient  number  of  representative  data  points  in  the  time  period  providing  the 
information  about  the  dynamics  (the  time  period  used  to  estimate  parameters).  Finally,  we 
note  that  the  model  developed  and  validated  here  retains  the  predictive  capabilities  present 
in  the  earlier  model  of  [3]  in  that  the  ability  to  predict  subsequent  data  after  fitting  the  model 
with  half  of  the  data  is  at  least  as  good  as  that  exhibited  in  [3]  for  the  various  patients  in 
our  data  sets. 


Patient:  Clinical  Indentification  Number  10 


Figure  12:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 
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In  Section  3  we  demonstrated  that  model  (2.1)  with  a  specific  example  set  of  parameters 
exhibited  multiple  equilibria.  Similarly,  after  fitting  the  full  set  of  patient  clinical  data  we 
used  the  parameter  sets  thus  obtained  to  determine  the  off-treatment  locally  a.s.  equilibria 
for  each  patient,  and  we  find  that,  for  all  but  two  of  the  patients,  there  is  only  one  lo¬ 
cally  a.s.  equilibrium  available  for  each  patient.  Obviously,  these  equilibria  are  all  infected 
equilibria,  with  equilibrium  viral  loads  ranging  from  450.2  to  51623  copies/ml-plasma.  In 
contrast,  patients  25  and  26  each  have  one  locally  a.s.  uninfected  equilibrium  and  one  locally 
a.s.  infected  equilibrium,  although  the  uninfected  equilibrium  for  both  these  cases  exhibit 
physiologically  unreasonable  levels  of  HIV-specific  CD8+  cells  as  compared  to  the  infected 
equilibrium  values.  This  behavior  is  not  unusual  since  we  do  not,  at  this  point,  have  any 
clinical  data  for  HIV-specific  CD8+  cell  counts  to  help  constrain  these  compartments.  In 
future  efforts  with  data  from  clinical  trials  currently  being  designed,  we  plan  to  explore  the 
use  of  an  additional  cost  term  in  the  optimization  algorithm  that  constrains  the  CD8+  values 
to  fit  observations. 


6  Conclusion  and  Remarks 


One  of  the  challenges  in  our  clinical  data  fitting  is  that  we  do  not  have  sufficiently  luxurious 
data  sets  to  fit  all  the  model  parameters  and  initial  conditions.  Due  to  the  high  correlations 
between  the  model  parameters,  biological  variations  in  all  parameters  across  population  and 
the  sensitivity  with  respect  to  parameters  varying  over  time  (especially  during  the  transition 
time  between  the  off-treatment  period  and  on-treatment  period),  we  intuitively  chose  some 
parameters  in  this  paper  to  be  fixed.  Even  though  the  simulation  results  indicate  reasonable 
fits  to  all  the  patients,  we  need  to  develop  a  scientific  methodology  to  deal  with  this  situation 
in  order  to  obtain  more  reliable  parameter  values  to  be  fixed  a  priori  and  thereby  sharpen  our 
estimation  results.  As  outlined  above  we  reduced  the  number  of  parameters  to  be  estimated 
by  fixing  some  of  the  model  parameters  at  population  averages.  In  clinical  cases  where  large, 
long  term  individual  patient  longitudinal  data  sets  are  not  available  as  in  our  study  here, 
one  could  use  population  averages  accrued  in  prior  patient  trials  to  reduce  the  number  of 
free  parameters.  Another  possible  method  to  reduce  the  dimension  of  parameter  space  and 
the  associated  high  correlation  between  these  parameters  is  principal  component  analysis, 
which  is  currently  under  investigation. 

The  clinical  data  fitting  results  in  Section  5  indicate  that  if  one  does  not  have  a  sufficient 
number  of  observations  during  the  periods  where  the  dynamics  are  changing,  then  it  is 
difficult  to  obtain  parameter  estimates  useful  in  predicting  corresponding  future  trends  in 
disease  progression.  Hence,  the  goodness  of  fit  results  are  affected  by  not  only  the  number 
of  observations  but  also  the  sampling  times  for  data  collection.  This  general  principle  was 
also  observed  and  explained  conceptually  in  [6]  where  the  general  “information  content”  in 
data  and  its  relation  to  sensitivity  is  explored.  These  considerations  suggest  that  parameter 
sensitivity  with  respect  to  data  is  important  in  experimental  design,  assisting  in  reducing 
effort  and  resources  required  to  collect  necessary  data.  One  potential  approach  for  this  is  the 
generalized  sensitivity  function  methodology  as  proposed  in  [36]  and  explored  in  [5,  7].  This 
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approach  combines  the  sensitivities  of  model  output  with  respect  to  model  parameters  with 
the  sensitivities  of  parameters  estimates  with  respect  to  changes  in  model  outputs.  Obtaining 
results  using  this  methodology  is  challenging  even  for  simple  examples  and  our  initial  efforts 
have  not  yet  proved  fruitful  for  our  complex  HIV  models  which  feature  high  correlations 
between  some  parameters  and  dynamics  that  dramatically  change  during  the  transition  time 
between  off-treatment  and  on-treatment.  Hence,  one  of  our  future  efforts  involves  develop¬ 
ment  of  a  methodology  to  determine  improved  sampling  times  for  data  collection.  We  are 
optimistic  that  ideas  contained  in  principal  component  analysis  and  generalized  sensitivity 
functions  can  be  combined  to  provide  new  guidance  in  this  regard. 
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Appendix:  Section  5  Results  for  Other  Patients 


Figure  13:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 


Patient:  Clinical  Indentification  Number  4 


Figure  14:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 
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Patient:  Clinical  Indentification  Number  5 


Figure  15:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 
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Figure  16:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 
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Patient:  Clinical  Indentification  Number  13 


Figure  17:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 


Patient:  Clinical  Indentification  Number  14 


Figure  18:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 


34 


Patient:  Clinical  Indentification  Number  24 


Figure  19:  Model  fits  to  data  (V)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 


Figure  20:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 
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Patient:  Clinical  Indentification  Number  26 


Figure  21:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 


Patient:  Clinical  Indentification  Number  27 
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Figure  22:  Model  fits  to  data  (V)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 
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Patient:  Clinical  Indentification  Number  33 


Figure  23:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 


Patient:  Clinical  Indentification  Number  46 


Figure  24:  Model  fits  to  data  (’x’)  with  parameters  estimated  from  half  longitudinal  data 
(solid  line)  or  full  data  (dash-dot  line).  Circles  denote  the  predicated  values  of  the  censored 
data  and  the  vertical  line  delineates  between  the  two  halves  of  the  longitudinal  data. 
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